Evidence for an oscillatory singularity in generic U(l) symmetric cosmologies on T 3 x R 
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A longstanding conjecture by Belinskii, Lifshitz, and Khalatnikov that the singularity in generic 
gravitational collapse is locally oscillatory is tested numerically in vacuum, (7(1) symmetric cosmo- 
logical spacetimes on T 3 x R. If the velocity term dominated (VTD) solution to Einstein's equations 
is substituted into the Hamiltonian for the full Einstein evolution equations, one term is found to 
grow exponentially. This generates a prediction that oscillatory behavior involving this term and 
another (which the VTD solution causes to decay exponentially) should be observed in the approach 
. to the singularity. Numerical simulations strongly support this prediction. 
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I. INTRODUCTION 

O ; 

An important open question in classical general relativity is the nature of the singularities that form in generic 
gravitational collapse. The singularity theorems of Penrose Q, Hawking ||||, and others Q|| prove that some type 
of singular behavior must arise in generic gravitational collapse of reasonable matter. However, these theorems do 
' not provide a description of the singular behavior that results. Many different types of singular behavior arise in 
known solutions to Einstein's equations. However, known explicit solutions tend to be characterized by simplifying 
symmetries. This means that singular behavior found in such examples need not be characteristic of those of generic 
collapse. In the 1960's, Belinskii, Khalatnikov, and Lifshitz (BKL) j6|-|Tot claimed to have shown that, in the approach 
to the singularity, each spatial point of a generic solution behaves as a separate vacuum, Bianchi IX (Mixmaster ( |Tl| ) 
homogeneous cosmology. Mixmaster cosmologies collapse as an infinite sequence of Bianchi I (Kasner p2]) spacetimes 
with a known relationship between one Kasner and the next [p|Jl^1 . Examinations of the BKL arguments (lj] and 
attempts to provide a more rigorous basis for their claims |Tq] have until recently yielded little evidence one way or 
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I , the other for their validity. 

■ In the past decade, work has begun aimed at understanding the singularity and the approach to it in spatially 
inhomogeneous cosmologies |lr^[l^1. If Einstein's equations are truncated by ignoring all terms with spatial derivatives 
and keeping all terms with time derivatives (for some more or less natural choice of spacetime slicing), the velocity 
k> term dominated (VTD) solutions are found. If an inhomogeneous spacetime is asymptotically VTD (AVTD), the 
evolution toward the singularity at (almost) every spatial point comes arbitrarily close to one of the VTD solutions 
[ 16, ^9|. It has been proven that polarized Gowdy cosmologies [ pfj|pl] ] and more general polarized T 2 -symmetric models 
[22] arc AVTD |l^,|2j|. The main extra feature in generic Gowdy models compared to polarized ones is the presence 
of two nonlinear terms in the Hamiltonian which yields the dynamical Einstein equations. In an AVTD spacetime, 
these must become exponentially small as the singularity is approached. But this requrement is only consistent with 
the form of the VTD solution if a spatially dependent parameter of the VTD solution lies in a restricted range at 
(almost) every spatial point. The numerical studies demonstrate that the nonlinear terms act as potentials to drive 
the parameter into the consistent range |p3. Very recently, it has been proven that generic Gowdy solutions with a 
consistent value of this parameter are AVTD [p5| . 

Recently, Weaver et al |2f| have extended the Gowdy model by inclusion of a magnetic field (and change of spatial 
topology from T 3 to the solv- twisted torus p7fl ). This model is an inhomogeneous generalization of the magnetic 
Bianchi VIo homogeneous cosmology which is known to exhibit Mixmaster behavior |28||29| ]. The magnetic field 
causes a third nonlinear term to be present. This term grows for precisely that range of VTD parameter which 
makes the original two vacuum Gowdy nonlinear terms exponentially small. This prediction of local Mixmaster 
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oscillations is observed in numerical simulations of the full Einstein equations. This study provided the first support 
in inhomogeneous cosmologies for the BKL claim. 

To further explore this issue, we have generalized to the study of vacuum spacetimes on T 3 x R with one spatial U(l) 
symmetry |30| . While such models are, of course, not generic solutions to Einstein's equations, they are considerably 
more complex than any previously considered for this purpose. Application of the methods previously described for 
Gowdy models leads to the predictions that (1) polarized U(l) symmetric cosmologies should be AVTD [^lj and (2) 
generic U(l) models should have an oscillatory singularity. Numerical simulations have previously provided strong 
support for the predicted AVTD singularity in polarized U(l) models |52|. Here we shall discuss the support we have 
obtained for the local oscillatory nature of the singularity in generic U(l) models. 

Unlike all previous cases discussed in this program ||l8| , p4 J26 32]| , numerical difficulties require an introduction 



of spatial averaging (data smoothing) at each time step to prevent numerical instability. This averaging destroys 
convergence of the solution with increasing spatial resolution (see |33|]). The need for spatial averaging can be traced 
to the growth of spiky features first seen and discussed in the Gowdy models |ll||24|]. These features also (as the 
bounces in Mixmaster itself Q]) make it necessary to explicitly enforce the Hamiltonian constraint. If this is not done, 
qualitatively incorrect behavior will result. We shall argue in this paper that, despite the need for spatial averaging 
and our simple algebraic method of enforcement of only the Hamiltonian constraint, the qualitative behavior of the 
numerical simulations is correct. We believe we have demonstrated that the singularity in generic U(l) symmetric 
cosmologies is spacelike, local, and oscillatory. Evidence for this is that the oscillations and their correlation with the 
exponential growth of certain nonlinear terms in Einstein's equations are independent of spatial resolution and choice 
of initial data. 

A brief review of generic U(l) models, their VTD solution, and the prediction of oscillatory behavior are given in 
Section II. In Section III we describe the numerical issues and results. Discussion is given in Section IV. 

II. THE MODEL 

As discussed in more detail elsewhere p2] |, U(l) symmetric cosmologies on T 3 x R are described by the metric 

ds 2 = e~ 2v {-e 2j V 4r dr 2 + e~ 2r e A e ab dx a dx b } + e 2ip (dx 3 + /3 a dx a ) 2 (1) 
where tp, j3 a , z, x, A are functions of spatial variables u, v and time r, sums are over a, b = u, v, and 



&ab — 9 



e 2z + e- 2z (l + x) 2 e 2z + e- 2z (x 2 -1) 
e 2z + e- 2z (x 2 -1) e 2z + e- 2z (l-x) 2 



(2) 



is the conformal metric of the u-v plane. A canonical transformation replaces the twists /3 a and their conjugate 
momenta e a with the twist potential u> and its conjugate momentum r |30 32 1. The dynamical variables ip and 



are respectively related to the amplitude for the + and x polarizations of gravitational waves and propagate in a 
background spacetime described by z, x, and A. Our coordinate choice (N = e A , zero shift) does not significantly 
restrict the generality of these models p0(| . Einstein's evolution equations (in vacuum) are found from the variation 
of |H 

H = I I dudvH 



dudv (\pl + \e iz p 2 x + \p 2 + ±e 4 V - \p\ + 2 PA ) 
+e~ 2T [ [ dudv{(e A e ab ) , ab -( e A e ab ) , a A, b +e A [(e~ 2z ) , u x, v - (e^ 2 ) , v x, u ' 



+2eVV,a <p,b +\e A e-^e ah u, a cu, b 



= Hk + Hy = J J dudvHx + J J dudvHy (3) 

where {p, r, p x , p z , p\} are respectively canonically conjugate to {ip, u>, x, z, A}. The constraints are 

H° = U - 2p A = (4) 

and 
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Ti-u = Pz z rU +p x x, u +p\ A,„ -Pa™ +P<P,u +ru m +5 { [e 4z - (1 + x) 2 ] p x - (1 + x)p z ) , v 

-±{[e 4z + (l-x 2 )]p x -xp z }, u =0, (5) 

H v = p z z, v +p x x, v +Pa A, v -Pa,v +P<Prv +ruj, v -\ { [e 4z - (1 - a;) 2 ] p x + (1 - x)p z ) , u 

+ \{[e iz + (l-x 2 )]p x -xp z }, v = {). (6) 

The VTD solution has been given elsewhere p2[ . Here we shall consider only the limit as r — > 00 of this solution. 
(Recall that the VTD solution is found by eliminating all terms with spatial derivatives from Einstein's equations.) 
The limiting VTD solution is 

z = -v z t, ,x = x , p z = -4w z , Px =p° x , ip = -v v t, 

v = u) 0) p= -4v v , r = r°, A = A + (2 - v a )t, p A = v A (7) 

where v z , v v , xo, p x , WOj ' , Ao, and va > are functions of u and v but independent of r. (The sign of va is fixed 
to ensure collapse.) 

We now use the method of consistent potentials (MCP) jlTj to determine the consistency of the VTD solution with 
the full Einstein equations. As r — > 00, consistency requires that all terms other than those which yield (Q) should 
become exponentially small. Rather than consider the equations, we shall examine the Hamiltonian density which 
generates them. Possible inconsistencies could arise from the nonlinear terms in (||) containing exponential factors. 
These same exponentials would, of course, be present in Einstein's equations. We first notice that the Gowdy-like 
terms fT|) 

V z = \ple^\ y^ir 2 e^ (8) 
in H.k become, in the limit of r — > 00, upon substitution of (R) 

V z ^\ V 2 x e-^>\ Vl ^\r 2 e- iv * T (9) 

and are exponentially small only if v z > and v v > 0. (As in the Gowdy case |24|| , non-generic behavior can arise at 
isolated spatial points where p x and/or r vanish.) 

The complicated terms in H containing the spatial derivatives have only two types of exponential behavior. All 
but one of the terms in 7iv have a factor 

e (-2T+A- 2z ) (1Q) 

(if we assume v z > 0, all components of e a b are dominated by e _2z ) which becomes 

m e (-*A+a«,)T (11) 

in the VTD limit. The remaining term is 

V 2 = 2 e e e w >a (12) 

which becomes upon substitution of (0) 

V 2 w F{x, Vo;)e ( -" A+2 "» +4 "» )T (13) 



where F is some function. The coefficients of r in ([Ll]) and (|l^) are restricted by the VTD form of the Hamiltonian 
constraint (as r — > 00) 

H° w -\v\ + 2v 2 z + 2v% ps (14) 

obtained by substitution of (0) into (Q). As discussed in (|lj) implies that va > 2w z so that ( [Tl| ) decays 

exponentially for v z > for any On the other hand, for V2 to become exponentially small with v z and v v > 0, we 
require v\ > (2v z + Av v ) 2 which is inconsistent with (p^). Since there is no way to make V\ and V2 both exponentially 
small with the same value of v v , the MCP predicts that either V\ or V2 will always grow exponentially. (Again, 
non-generic behavior can result at isolated spatial points where the coefficient of V2 happens to vanish.) 
To refine this prediction, consider v z » \v v \. Substitution in ( p^ ) yields 
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v A ^2v z + ^ (15) 

which shows that 

V 2 « e 4v « T F(x, Vw) , Viwr 2 e" 4 ^ T . (16) 

Thus we expect V\ and Va to act as potentials for the (p degree of freedom with a bounce off cither potential causing 
the sign of v v to change. The remaining variables will follow the VTD solution with parameters which change at 
every bounce in (p. Thus we predict that oscillations in the ip degree of freedom will occur at (almost) every spatial 
point with different values of v v and coefficients of V\ and V2 ■ 

Note that in polarized U(l) models, u> = r = so that the oscillations of ( |l6| ) should be absent. Polarized U(l) 
models should thus be AVTD. This is precisely what has been found in numerical simulations of these models p2fl. 



III. RESULTS AND NUMERICAL ISSUES 



In order to test the predictions of the previous section, we performed numerical simulations of the Einstein evolution 
equations obtained from the variation of ||). We use a symplectic PDE solver which has been described in great 
detail elsewhere |l^,|3^,[55| ]39| . In our previous study of polarized U(l) models, we demonstrated convergence of the 
solutions at the expected order with increasing spatial resolution. Unfortunately, the re-introduction of the u> degree 
of freedom leads to the growth of spiky features absent in the polarized case. The origin of the spiky features is 
discussed elsewhere [|4| but is related to the non-generic behavior at isolated spatial points. The spiky features which 
our methods easily treat in one spatial dimension []l8f cannot be modeled sufficiently well in two spatial dimensions 
to prevent numerical instability. However, these instabilities can be suppressed in two ways. First, the simulations 
are begun at r « 10 or so to reduce the influence of TLy from (H|). (This has the disadvantage, as in homogeneous 
Mixmaster spacetimes p0[ , of increasing the time interval between bounces.) Second, at every time step, all ten 
variables are replaced by their spatial "average." Any function f(u, v) is replaced by ||3S| ] 
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f(u, v) = C f(u, v) + ^ C i if( u + v ) + f( u > v + i^ v ) 

i=l 

+f(u - iAu, v) + f(u, v - iAv) - 4/(u, v)} (17) 

where c = 1/2, c x = 4867/38400, c 2 = -1067/28800, c 3 = -1237/691200, c 4 = 787/345600, c 5 = 31/691200 and 
Au, Av are the grid spacings. This scheme is 6th order accurate (i.e. the difference between f(u,v) and f{u,v) is 
7th order in the grid spacings — actually 8th order due to symmetry). Where f(u,v) is smooth, the averaging has no 
effect. However, it does spread out grid scale size spiky features. Unfortunately, the averaging process destroys the 
convergence seen in the polarized case (by examination of the deviation of the Hamiltonian and momentum constraints 
from zero j33|), as can be determined by adding averaging to polarized C/(l) simulations. We shall argue later that 
despite the absence of convergence in this sense, the qualitative behavior is still independent of spatial resolution. 
We shall further argue that the correct qualitative behavior is sufhcent to determine the nature of the singularity in 
generic U(l) models. 

In Section II, we showed that the behavior of the potentials in TLy was restricted by the VTD limit of the Hamiltonian 
constraint ( pT|) . This means that it is essential to preserve the Hamiltonian constraint during the simulation. Failure 
to preserve 7{P = means that p\ would have the wrong dependence on p z and p. This could change the sign of the 
coefficient of r in (|l5|). While solving the constraints initially is sufficient analytically, there is no way to guarantee 
that the differenced form of the constraints is preserved during a numerical simulation of the Einstein evolution 
equations. While in the polarized U(l) case, the constraints can be seen to converge to zero p2] |, this is no longer 
true in generic U(l) simulations with spatial averaging. To preserve the correct relationships among v\, v z , and v v 
from (|l4|) , we therefore solve TL° = (Q) for p\ > algebraically at each time step. We note that this is the precise 
analog of the procedure used in Mixmaster itself [j34| . However, spatial averaging must be performed on this new p\ 
which then yields a (small) non-zero value for Ti. . If we assume that pa = p\ IUC + A is the measured value, then ( |l4| ) 
can be linearized in the error A and solved to yield 

A = -LL (18) 

PA 

where H. is the measured, erroneous value of the Hamiltonian constraint. Substitution in a = —v\ + 2v z + 4w v , the 
coefficient of r in V2, yields a measured and true value for a. If the product aa true is positive everywhere for all r, 
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the errors due to averaging cannot have changed the qualitative behavior since that just depends on the sign of a. 
Examination of the simulation data shows a a true > at (almost) all spacetime points. Occasionally, early in the 
simulation, a few isolated points — in space and time — will show a sign change due to this error. 

The momentum constraints (^) and (0) are freely evolved and do not remain especially small. However, they 
contain only spatial derivatives. Errors in the momentum constraints would therefore generate errors in the spatial 
dependence of the variables at a given time but not in their qualitative time dependence at each spatial point. A 
measure of the constraint convergence vs r is shown in Figure 1. The freely evolved TL U (Hv ~ Ti. u ) is actually larger 
at finer spatial resolution. This may be attributed to the larger spatial gradients observed at finer spatial resolution 
. On the other hand, the error in the Hamiltonian constraint is converging to zero. 
The restricted initial value solution for generic U(l) models has been described elsewhere p2|| . To solve the 
momentum constraints, we set p z = p x = <p, a = uj ia — and p\ — ce A . For c > sufficiently large, this allows the 
Hamiltonian constraint to be solved algebraically for p (case A) or r (case B). Four functions, x, z, A, and r (case A) 
or p (case B), may be freely specified. The detailed values used in the numerical simulations are given in the figure 
captions. While we shall present results for only two sets of initial data, the results are typical of all initial data and 
may be regarded to be characteristic of generic initial data. 

Figures 2-4 show the behavior of ip, Vi, and V2 for three spatial resolutions for case A initial data. In all cases, the 
behavior is exactly as predicted in Section II. Figure 5 shows a typical evolution for case B initial data, again showing 
the predicted behavior. For the simulation of Figure 3, the remaining variables lo, z, x, and A vs r are shown at the 
same spatial point in Figure 6. Here we see approximate VTD behavior inw, z, and A with parameters changing at 
the times of the ip bounces. While the behavior of x does not appear to be as predicted, we note from (Q), (|l0|), and 
( |l2| ) that only the behaviors of z and A are important to the local dynamics since only they appear in the arguments 
of exponentials. Finally, we display "movie frames" of TCy(u,v) and ip(u,v) vs r in Figures 7 and 8. S ince 1~Ly is 
shown on a logarithmic scale and tp on a linear one, Figure 7 is dominated by V% while Figure 8 indicates the behavior 
on a logarithmic scale of V\. The predicted behavior that either V\ or V2 but not both be large is clearly seen in 
Figures 7 and 8. The bounces in (p occur at different spatial points at different times. This eventually will lead to 
increasingly complex spatial structure E]] . 



IV. DISCUSSION 



Numerical simulations of generic U(l) models demonstrate that the evolution toward the singularity is local and 
oscillatory due to the alternate expoential growth and decay of V± and V2. This behavior requires that the time 
dependence of the coefficients of V\ and V2 and the coefficients of r in the approximate forms of V\ and V% be 
negligible (i.e. act on a much slower time scale) compared to the exponential time dependence obtained from the 
VTD solution. The qualitative behavior seen in the simulations indicates that this requirement is met (almost) 
everywhere sufficiently close to the singularity. Exceptional behavior at isolated spatial points similar to that found 
in simpler models p4| , p6[ cannot be studied with the current numerical code. Exceptional behavior arises at isolated 
spatial points because one of the exponential potentials which causes the generic behavior is absent. While we believe 
we have identified the correct exponential behavior, the data we have on spatial dependence at a given time is probably 
not sufficiently reliable for conclusions about "higher order" effects. 

The explanation for the local nature of the evolution is easy to obtain. From the metric (|l|), we see that the distance 
AZ traveled by a light ray away from the singularity from r = 00 to r = tq (coordinate horizon size) is 

/•TO 

Aira / e A/2+z ~ T d T (19) 
J 00 

where the VTD limit of e a j, has been used. Simulations demonstrate that the VTD solution may be used in ( |l9| ) to 
give AZ — > for the coordinate horizon size pl[| as To — ► 00. Since the horizon size is decreasing, the spatial points 
are unable to communicate with each other. But such communication occurs through changing spatial derivatives. If 
no communication can occur, the spatial derivative containing terms must be dynamically unimportant. 

The primary question then is the extent to which the numerical results presented here are believable as evidence 
for the oscillatory nature of the singularity. It is easy to show, e.g. by comparing polarized J7(l) simulations with 
and without spatial averaging, that spatial averaging ruins convergence tests. Furthermore, non-zero values of the 
momentum constraints (^) and (^) indicate the presence of errors in the variables and their spatial derivatives. (Since 
rescaling the coordinates by a constant yields rescaled values of the constraints without changing the behavior of the 
solutions, one cannot attach any significance to the actual magnitude of the constraint violation.) Both sources of 
error — spatial averaging and constraint violation — principally affect details of the spatial dependence of the variables 
at a given time. Comparison of simulations at different spatial resolutions (and with quite different initial data) shows 
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that the qualitative time dependence at fixed spatial points is not affected by these errors. This is shown in Figure 9 
where tp is shown on a line u = v with u, v G [0, tt] vs t for three different spatial resolutions. While the features appear 
narrower at finer spatial resolution (see p3 for a similar phenomenon in Gowdy models) , the time development of the 
simulations is remarkably consistent and independent of resolution. Figure f shows the final time step for tp for two 
different spatial resolutions. While the spatial dependence is quantitatively different, it is qualitatively quite similar. 
We also see in Figures f f and 12, that the correlation between V\ and Vi seen at single spatial points in Figures 2-6 
occurs everywhere. Figure ff shows log 10 V\ and log 10 V% on the line u — v with u,v € [0, 2ir] vs r. Clearly, one 
potential grows as the other decays at all spatial points. Figure 12 shows the final time step for Vt, V2, and tp on the 
line u = v with u, v S [0, n]. Given this same oscillatory behavior in all the simulations and the explanation for it 
in terms of the VTD solution given in Section 11, it is proabable that we have correctly described the nature of the 
generic singularity in this class of models. 

This does not mean that it is unimportant to improve the numerical treatment. Better numerics will yield more 
accurate spatial dependence of the metric variables. We expect to see ever narrowing spatial structure as in Gowdy 
and magnetic Gowdy [^6| models caused by non-generic values at isolated spatial points. This can also yield 
precise data for analysis of the local dynamics to characterize the "higher order" effects due to spatial inhomogeneity. 
ft might also become possible to observe and characterize any other exceptional behavior that might arise at isolated 
spatial points. 

There are several possibilities for numerical improvement. First among these is to incorporate into the algorithm 
the fact that a known explicit solution exists for the bounce (scattering) off an exponential potential. The symplectic 
algorithm divides the Hamiltonian for a system into two or more subhamiltonians with explicit exact solutions for 
their equations of motion. The standard division of H is into kinetic (H^) and potential {Hy) pieces. Including a 
dominant exponential wall from Hy in Hk has yielded tens of orders of magnitude of improvement in speed while 
maintaining machine level precision in the ODE'S of spatially homogeneous Mixmaster models [ [34| . The current U(l) 
code uses TLk and TLy (|^) as the subhamiltonians. However, one could treat 

= \p 2 + \e K e ab e^u, a u, b (20) 

or 

n^ = \ P 2 + ye^ (21) 

as a separate subhamiltonian depending on which of V\ , V2 is the larger. The absence of pa, Pz, and r from (|2^) means 
that the equations from Tip are exactly solvable (as, obviously are those from Ti^). While in the ODE case, the 
main advantage was found in the ability to take huge time steps, for the U(l) case, it will be the improved accuracy 
of the bounce solution. Another area of improvement is in spatial differencing. The current code uses 4th and 6th 
order accurate representations of first and second derivatives due to Norton |3^] . These are designed to minimize the 
growth of grid scale wavelength instabilities. However, when applied to the Gowdy model, this scheme appears to 
be less stable than the original one Jl8| based on variation of a differenced form of the Hamiltonian. We also note 
that spectral methods have been tried but have not proven to be useful. A third place for numerical improvement 
would be to solve all three constraints rather than just the Hamiltonian constraint. Of course, one naturally wonders 
if adaptive mesh refinement (AMR) might yield improvements in accuracy and stability. Studies to date have shown 
that increasing spatial resolution yields only small improvements in stability. This is probably due to the fact that 
finer resolution only gives better representation of spiky features — showing them to be narrower and steeper than 
they appear at coarser resolution. In the Gowdy simulations, however, it was noticed that, for any given simulation, 
there is a threshold spatial resolution. If the resolution is coarser than that, the code will blow up before the AVTD 
regime is reached everywhere. It is therefore possible that we have not yet reached this threshold in the U(l) case 
and thus would be helped by AMR. However, in the generic U(l) models, there is no AVTD regime. Thus it is not 
clear if any spatial resolution could yield a simulation which could run indefinitely without crashing. 

As a final note, we remark that the local nature of the U(l) evolution (as well as the other cases we have studied) 
has two major implications: 

(1) The use of cosmological boundary conditions is merely a convenience since they do not affect the local behavior. 
Any collapse of a system with one spatial Killing field to a spacelike singularity should be local and oscillatory. 

(2) Qualitative answers on the nature of the singularity do not require fine spatial resolutions. This means that the 
zero Killing field case should be tractable numerically. Work on this case is in progress. 



6 



ACKNOWLEDGEMENTS 

We would like to thank the Albert Einstein Institute at Potsdam for hospitality. BKB would also like to thank the 
Institute for Geophysics and Planetary Physics of Lawrence Livermore National Laboratory for hospitality. This work 
was supported in part by National Science Foundation Grants PHY9507313 and PHY9503133. Numerical simulations 
were performed at the National Center for Supercomputing Applications (University of Illinois). 



7 



FIGURE CAPTIONS 



Figure 1. Convergence testing of constraints. The average values of the momentum constraint Ji u (broken line) 
and Hamiltonian constraint H° (solid line) are displayed vs r for generic vacuum U(l) symmetric model simulations 
with 256 2 (nothing) and 512 2 (triangles) spatial grid points. (See Figure 2 for initial data.) 

Figure 2. Oscillatory behavior at a typical spatial grid point. The potentials V\ and V2 and ip are shown vs r for 
a simulation with 128 2 spatial grid points. The initial data are A = sinu sinv, x = z = cosu cosw, ip = uo = 0, 
Pa = 14 e A , r — 10 cos u cosv, and p z = p x = 0. The Hamiltonian constraint is solved for p. The initial value of t is 
10. 

Figure 3. Oscillatory behavior at a typical spatial grid point. The same initial data as in Figure 2 but with 256 2 
spatial grid points. 

Figure 4. Oscillatory behavior at a typical spatial grid point. The same initial data as in Figure 2 but with 512 2 
spatial grid points. 

Figure 5. Oscillatory behavior at a typical spatial grid point. The same as Figure 2 but with Case B initial data: 
A = sinu sin?;, x = z = cosu cosv, p = uo = 0, pa = 4e A , p = cosu cosv, and p z — p x = 0. The Hamiltonian 
constraint is solved for r. The initial value of r is 10. 

Figure 6. Behavior of (a) A and z, (b) u, and (c) x in the simulation of Figure 3. 

Figure 7. Movie frames of log 10 Hv (u, v) at (from right to left and top to bottom) r = 14.9, 19.8, 24.7, 29.6, 34.5, 
39.5, 44.4, 49.3, 54.2. The range of values is —30 (black) to 3 (white). The simulation of Figure 3 provided the data. 

Figure 8. Movie frames of ip(u,v) at the same values of r as in Figure 7, taken from the same simulation. The 
range of value is —27 (black) to 3.8 (white) on a linear scale. 

Figure 9. The influence of spatial resolution. The variable <p on the line u = v for u,v G [0,7r] (vertical axis) is 
shown for r G [10,53.5] (horizontal axis). The simulations of Figures 2 (top), 3, and 4 (bottom) are shown. The gray 
scale is similar to that in Figure 8. 

Figure 10. The variable ip is shown for the line u = v for u,v G [0,7r] at t = 53.5 for a simulation with 256 2 
spatial grid points (broken line) and 512 2 spatial grid points. These are the same simulations as in Figures 3 and 4 
respectively. 

Figure 11. LogioFi (left) and log 10 V2 (right) for the line u = v for u, v G [0, 27r] (vertical axis) and r G [10,59] 
(horizontal axis) are shown on the same scale. The simulation of Figure 3 is used. 

Figure 12. The final time step of Figure 11 is shown for log 10 V\ (broken line), log 10 V2 (thick solid line), and <p 
(thin solid line). The horizontal axis is the line u = v for u, v G [0, n] 
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